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Abstract 

Filon-Simpson quadrature rules are derived for integrals of the type 

Jadx f{x) sin{xy)/{xy) and jj^ dx f (x) A sin^ {xy /2) / (xy)"^ 

which are needed in applications of the worldline variational approach to Quantum Field Theory. 
These new integration rules reduce to the standard Simpson rule for y = and are exact for y —>■ oo 
when a — and /(O) ^ 0. The subleading term in the asymptotic expansion is also reproduced 
more and more precisely when the number of integration points is increased. Tests show that the 
numerical results are indeed stable over a wide range of j/-values whereas usual Gauss-Legendre 
quadrature rules are more precise at low y but fail completely for large values of y. The associated 
Filon-Simpson weights are given in terms of sine and cosine integrals and have to be evaluated for 
each value of y. A Fortran program to calculate them in a fast and accurate manner is available. 



1 Introduction 



^0 



Many problems in theoretical and computational physics require numerical integration over rapidly 
oscillating functions which usually is considered a difficult, if not hopeless problem. However, as 
Iserles has pointed out in ref. [1] this must not be the case: If the quadrature rule for an integral of 
the type 

rh N 

dx/(x)e^"^ ~ (1.1) 

1=0 

includes the endpoints xq = 0,xn = h then the error can be reduced to O {h^^'^ /y^^ for y — oo. 
In Iserles' words: As long as right methods are used, quadrature of highly-oscillatory integrals is very 
accurate and affordable ! This result rules out the usual quadrature schemes of Gaussian type but is 
realized in the classic Filon integration rules [2] which are covered in many standard books (see, e.g. 
chapter 2.10 in ref. [3] or eqs. 25.4.47 - 25.4.57 in ref. [4]). 

The present note is not concerned with the numerical evaluation of the Fourier integral (1.1) which 
is well treated in mimcrous works but with oscillatory integrals whose integrand has a removable 
singularity at the origin, viz. 

h[f]{a,b,y) := f'dxfix)^-^^ (1.2) 
J a xy 



h[f]{aAy) := t dxf{x)A^-^^^yP-. (1.3) 
Ja X y 



and 

sin^(x|//2) 

x'^y'^ 

The corresponding weight functions have been normalized in such a way that they become unity for 
vanishing external (frequency) parameter y. In eq. (1.3) one could rescale y ^ 2y to obtain a simpler 
oscillating weight function. However, due to sin^ t = (1 — cos(2t)/2 it would oscillate with twice the 
frequency and therefore we prefer to define I2 in the above form. This is also how it appears in the 
applications we will discuss below. Note the relation 

7i[/](a,6,y) = [y^hifKaAv)] ■ (1-4) 

As frequently is the case the necessity to evaluate numerically integrals of the type (1.2, 1.3) arose 
from a concrete theoretical problem. Here it was the polaron problem in solid state physics [5] and 
its extension, the worldline variational approach to relativistic Quantum Field Theory [6, 7, 8, 9, 10, 
11, 12]. With a general quadratic trial action the Fcynman- Jensen variational principle gives rise to a 
pair of non- linear variational equations for the "profile function" A[E) and the "pseudotime" //^(cr). 
For example, when calculating the self-energy of a single scalar particle in d space-time dimensions 
these equations are of the form 

1 

ME) 



2^ ^ 4 /^^p 1 sin"(£;a/2) 



2 



. (1-8) 



Here V[iJ,^] is the interaction term averaged over the trial action. It is a functional of the pseudotime 
and specific for the field theory under consideration. For small proper times a its functional derivative 
(or "force") has the behaviour 

6V CT^o constant 

where r = for super-renormalizable theories and r = 1 for renormalizable ones like Quantum Electro- 
Dynamics (QED). Thus for the 3-dimensional polaron problem and for super-renormalizable theories 
in 4 dimensions no divergences are encountered in the variational equations but QED4 needs extra 
regularization, e.g. a cut-off at high momentum/small proper time. In any case, at least the relation 
(1.6) between profile function and pseudotime requires evaluation of integrals of the type (1.3) and 
frequently the proper-time derivative 

da ~ TT io A{E) 'W 
i.e. integrals of the type (1.2) are also needed. 

Of course, eqs. (1.5, 1.6, 1.8) involve infinite upper limits but the large- i?, cr limit of profile 
function and pseudotime may be obtained analytically so that we can restrict ourselves to finite 
integration limits and add the analytically calculated asymptotic contribution to the numerical result. 
Numerical evaluation of Fourier integrals with infinite upper limits is much more demanding; one 
strategy (employed in the adaptive NAG routine DOIASF [13] based on ref. [14]) requires a delicate 
extrapolation procedure and thus seems to be not suitable for a numerical solution of the variational 
equations. 

Previously we have solved these equations on a grid of Gaussian points by iteration (for details 
see refs. [7]) using Gauss-Legendre quadrature for the oscillatory integrals. Although sufficient for 
many purposes problems of numerical stability became more serious, in particular in QED4 when the 
momentum cut-off was increased too much. It should be mentioned that in the 1-body sector the 
profile function can be eliminated alltogether giving rise to a non-linear integro-differential equation 
for the pseudotime only which shows a striking resemblance to the classical Abraham-Lorentz-Dirac 
equation [15]. This eliminates the need for numerical evaluation of oscillatory integrals but requires 
solution of non-linear delay-type equations. In addition, at present this is restricted to the 1-body 
case and for the 2-body, bound-state case we had to apply the previous scheme based on (several) 
profile functions and pseudotimes [12, 16]. Hence there is a definite need to have a reliable, fast and 
stable method for evaluating integrals of the type (1.2, 1.3). 

In the following sections these integration rules are derived "naively" (that is without mathematical 
rigor and error estimates) and tested for simple cases. Of course, there is a rich mathematical literature 
on quadrature rules for various or general oscillatory integrals (for recent developments see, e.g., refs. 
[17, 18, 19]) but here the emphasis is more on the practical implementation and availability. In 
addition, the behaviour for large values of the parameter y will be investigated in some detail. 

2 Filon- Simpson quadrature 

The strategy for deriving a stable quadrature formula for oscillating integrals of the type 

l-l, N ( sin{xy) . j = I 

/ dxf{x)Oj{xy)c.J24'^f{4'^) , 0,{xy) = 4sin^(.7/2) (2-1) 
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is well known [1]: choose N points in the interval [a, 6], two of them identical with the end- 
points and require that the integral over x^Oj{xy) is exact. This gives a system of equations for the 
integration points and weights wl '. 



Here, for simplicity, we choose equidistant, j-independent points 



x^P = Xi = a + ih , i = . . . N , h = —^y— 



(2.2) 



and N = 2. In other words: this will be a generalization of Simpson's time-honoured rule, to which it 
should reduce for Oj{xy = 0) = 1. It is for this reason that we may call it Filon-Simpson quadrature. 
The integration points being fixed there are three monomials which can be integrated exactly giving 
rise to three equations for the weights 



=a+2/i 



Y^w^px\ = [ '^'^'^dxx'^Ojixy) := k = 0,l,2 



It is easy to solve this linear system of equations with the result 

XI X2 J^^ - (xi + X2) Jp'^ + J^^ 



w 





1 r 




2h? - 


? = 


1 r 






= 


1 r 




2^ 



-2x0 X2 Jq'^ +4x1 Ji ^ - 24^^ 



Xq xi - (xo + Xi) J}^^ + 



rU) 



(2.3) 

(2.4) 
(2.5) 
(2.6) 



Note that the weights depend on the lower and upper integration limit as well as on the external 
parameter y: 



= w'f\a,b,y) , h 



The relation (1.4) translates into 



J, 



(1) 



1 d 



2y dy 

and the integrals may be written as 



2/2 jf, fc = 0,1,2 



w. 



(1) 



J, 



with 



k+l 



Fl'\hy) - Fi^\ay) 



dtt^Oj{t) . 



b — a 



(i) 

Prom the explicit form of the weight functions Oj: one finds 



F^'^\z) = Si{z) , Fl"'{z) = 1-cosz, F^"'{z) = sinz-zcosz 



.(1) 



1 — cos z 



, Fi'^\z) = 2 [7 + lnz-Ci(z)] 



2 [z — sin z ] . 



(2.7) 
(2.8) 

(2.9) 
(2.10) 

(2.11) 
(2.12) 
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Here 7 = 0.5772156640 ... is Euler's constant and Si{z),Ci{z) the sine and cosine integral, respectively, 
as defined in chapter 5 of ref. [4]. Note that the integrals jj.j'^ are odd under exchange a <-> 6 as may 
be seen from their definition or from eq. (2.9). Therefore we find 

w\'\b,a,y) = -w'i\{a,b,y) , i = 0,1,2 (2.13) 

which reflects the basic property of the integrals Ij[f] when their limits are exchanged. 

As a check we now consider the limit y ^ 0. It is easily seen that the integrals jj^^^ and therefore 
the weights are well-behaved in that limit because we have 

00 • 2n 

F»w = .-|:,-ir^^-^. (2.14) 



Therefore 



'^'^ k + l 



Eqs. (2.6 - 2.4) and b = a + 2h then immediately give the standard Simpson weights 

+ (2.16) 

independent of the lower limit a and the type j of the integral as it should be. 

More interesting is the limit y ^ 00 for the case a = 0. First, the exact asymptotic behaviour of 
the integrals Ij is easily obtained by substituting t = xy : 

/,[/](« = 0, 6, y) = i ["'dtfC-) 0,it) M rdtO,{t) = M^, ,- = 1,2. (2.17) 

y Jo \yJ y JO y 2 

The asymptotic limit of the weights is found by employing the expansions [4] 



which leads to 



Q., X z-*oc VT cosz smz z-^oo smz 

Si(^) ^ 2 ~ ~z + • • • ' Ci(z) ^ + ■ ■ ■ (2.18) 



J? %^l-^ + o(i^| (2.19) 



Therefore Jq^ provides the leading asymptotic contribution and from eqs. (2.4 - 2.6) we see that for 



xo = a = only Wq' is affected by it: 



0) j ^ (2i -2) \nby,l,cosby 

y2+^V y 



^^^^ y.j-.,...uy,.,.^.uy ^ (2.20) 



0(^(2i-2)ln^^y,l,cos6y-^ ^2.21) 
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Hence Filon-Simpson quadrature gives in that limit 

jZwi^^fiih) ^-^4^V(0) = ^f/(0) (2.22) 

in perfect agreement with the leading asymptotic result (2.17) ! 

Of course, the story is different if the function f{x) vanishes at x = 0: then the next-to- leading 
terms of the asymptotic expansion come into play. These are studied in more detail in Appendix Al 
and A2 where it is shown that for small enough increment h the subleading terms are also reproduced 
approximately. For example, derivatives are replaced by finite differences as is expected from a 
quadrature rule which is based on values and not derivatives [19] of the function to be integrated. 



3 Extended Filon-Simpson rule 

We would like to decrease the error of the quadrature rule in a systematic way without messing up 
the construction. For equidistant integration points this is, of course, very easy to achieve by dividing 
the integration interval in N jl sub-intervals, being even, and applying the integration rule in each 
sub-interval. We then obtain the 'extended' (or 'composite') Filon-Simpson quadrature rule 

/ dxf{x)Oj{xy) ~ ^ VF^^^V(« + , h = — ^ , A even (3.1) 

where (using the explicit nomenclature of cq. (2.7)) 

W!^'^ = wi^\a,a + 2h,y) (3.2) 
Wp^ = w[^\a + {i-l)h,a + {i + l)h,y) i = 1, 3 . . . (A - 1) (3.3) 
Wj;^^ = w^^\a+{i-2)h,a + ih,y) + wi^^\a + ih,a + {i + 2)h,y) i = 2, 4 . . . (A - 2) (3.4) 
Wij^ = wi^\a+{N -2)h,a + Nh,y) . (3.5) 

(i) 

Because of the complicated dependence of the weights w^^ on the lower and upper integration limit 
it is not possible to simplify the above expressions for y 0. Only for y = we obtain the weights 
for the extended Simpson rule 

: i = 0,N 

: f = l,3...(A-l) (3.6) 
: i = 2,4... (AT -2). 

The asymptotic limit y ^ oo for a = is also unchanged from the the previous section: this is because 
only in the first interval from a = up to 2h the leading term 7r/2 from the asymptotic expansion 
(2.18) of the sine integral survives whereas in the next and all other intervals it is cancelled due to 
difference which has to be taken in eq. (2.9). Thus in the extended Filon-Simpson quadrature rule 
asymptotic high oscillations are also treated correctly in leading order. Appendix A3 shows that this 
also holds for the subleading term if the increment h is made small enough, i.e. if the number N of 
integration points is increased. This is a remarkable feature since logarithmic terms show up in the 
next-to-leading order of the asymptotic expansion of l2[f] (see Appendix Al). However, eq. (2.12) 
shows that by construction such terms are also present in the Filon-Simpson rule for I2 [/] . 
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4 Numerical tests 



Here we report test evaluations of the integrals 



I, 



fi = x^e (^a = 0,b = oo,y^ = f dx e Oj{xy) 



(4.1) 



for / = 0, 1 with the Filon-Simpson and other quadrature rules. The oscillatory weight functions 
Oj{xy) have been defined in eq. (2.1). The exact values of the test integrals are (see e.g. p. 234 - 
235 in ref. [20]) 



h[k] 



arctan(7/) 



y 



h[fi 



1 



2yarctan(2/) -ln(^l + y^jj , hifi] = ^ln{l + y' 



(4.2) 
(4.3) 



and fulfill the relation (1.4). As discussed in the Introduction our strategy is to perform the numerical 




Figure 1: Absolute value of the relative deviation between numerical and exact result for the oscil- 
latory integral /i[/o] as a function of the frequency parameter y. The numerical result was obtained 
with different quadrature rules by numerical integration up to 6 = 20 and adding the asymptotic 
contribution (4.6). The integrand is the function fo{x) = times the oscillatory weight function 
Oi (xy) = sin(xy) / (xy) and 144 integration points have been used. 
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Figure 2: As in Fig. 1 but with 288 integration points. Note the extended scale for the relative 
deviations. 



integration up to x = 6 and to add the asymptotic contribution 



dxx^ e ^ Oj (xy) . 



(4.4) 



Due to the specific form (2.1) of the weight functions Oj it can be expressed in terms of the exponential 
integral [4], viz. 



6' 

Ai, = -lmEi_i{{l-iy)b) , A21 = 

y 

The asymptotic expansion of these functions gives 



26'- 



Re 



E2_Kfo) - E2_/ ((1 - %)6) 



^21 - 2- 



l + y2 



sin by / 1 

cos fey + - + 0(- 

y \b 

1 + y"^ + y sin by — cos by / 1 



+ 



T + y2 

which is well-behaved even for y = 0. In all numerical calculations we have chosen 

b = 20 



(4.5) 

(4.6) 
(4.7) 

(4.8) 
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so that in general only a very small fraction O (e^^*') ~ 5 • 10^^ times powers of b comes from the 
asymptotic region. We have checked that the next order of the asymptotic expansion does not alter 
the outcome on the required accuracy level except at y = where small changes are observed. 

The Filon-Simpson weights are given in terms of sine and cosine integrals so that a fast and 
precise routine for these special functions is mandatory. Fortunately there are many convenient 
rational approximations available, e.g. in chapter 5 of ref. [4] or routines from commercial libraries. 
However, with little extra effort the Chebyshev expansions given in Table 23 of Luke's comprehensive 
book [21] provide accurate values for these functions which exceed the double precision arithmetic 
which we employ. Therefore we have used these expansions in a Fortran subroutine making the whole 
program for calculating Filon-Simpson weights self-contained and portable ^. 

Fig. 1 shows the relative deviation of the numerically calculated integral /i[/o] plus asymptotic 
contribution from the exact result (4.2) as a function of the external variable y. 144 integration 
points have been used for both the Filon-Simpson and Simpson rule and 2 x 72 points for the Gauss- 
Legendre quadrature (i.e. subdivision of the whole interval in 2 parts and application of a 72-point 
Gauss-Legendre rule in each part). In Fig. 2 the number of integration points is increased to 288 
and 4 X 72 respectively. It is seen that for small y the Gauss-Legendre rule is superior but starts 
to deteriorate when the number of integration points is insufficient for the rapid oscillations of the 
integrand. 
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Figure 3: Same as in Fig. 1 but for the test function fi{x) = xe ^ and 288 integration points. 



copy of the program together with a sample run is available on request. 
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By increasing the number of integration points the onset of failure can be extended (from y ~ 25 
in Fig. 1 to y ~ 50 in Fig. 2) but not avoided ^. The ordinary Simpson rule is even less capable 
to deal with such type of integrals. In contrast, the new Filon-Simpson quadrature rule gives stable 
results for all y- values considered and its relative deviation from the exact result is a smooth function 
of y decreasing for very large y. Of course, the calculation of the weights has to be redone for each 
value of y and is more involved than for the standard, simple integration rules. However, in terms 
of CPU-time this is still a neglible expense (200 y-values in Fig. 2 took about 2.5 seconds on a 600 
MHz Alpha workstation). Fig. 3 depicts the result for the test function /i which has an additional 
x-power in the integrand so that the relative accuracy which can be achieved with a fixed number of 
integration points is worse than in the previous case. In addition now /i(0) = so that the leading 
asymptotic term vanishes but Appendix A3 demonstrates that the subleading term is also reproduced 
for sufficiently large number of integration points. Therefore Filon-Simpson integration still does far 
better than ordinary Simpson or Gauss-Legendre rules. The corresponding results for the integral 
/2[/z], i.e. with the weight function As\v?{xy /2) / {xyY are shown in Figs. 4 and 5 and confirm the 
experience gained with Ii. 




100 



Figure 4: Same as in Fig. 2 but with the oscillatory weight function 02{xy) and the asymptotic 
contribution (4.7). 



rule-of-thumb is that one needs at least one Gaussian integration point on each oscillation. Indeed, assuming 
that values up to a; = 6 contribute significantly to the integrals we would have ymax — nN/b ~ 0.15A'' in qualitative 
agreement with Figs. 1, 2. 
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Figure 5: Same as in Fig. 3 but with the oscihatory weight function 02{xy) and the asymptotic 
contribution (4.7). 



Finally, we have investigated whether it is advantagous to use the relation 

h[f]{aAy) = \ r dy'y'h[f]{a,b,y') (4.9) 

y Jo 

which is obtained from eq. (1.4) by integration (the integration constant is zero since the integrals Ij 
are finite at y = 0). Here one doesn't have to integrate over an oscillating function and the asymptotic 
limit (2.17) is also correctly obtained. In the worldline application this would amount to first evaluate 
dfi'^{a)/da from eq. (1.8) and then integrate it step by step via a trapezoidal or Simpson rule to 
obtain fj,'^{a). 

Fig. 6 shows the comparison with the direct Filon-Simpson integration. While a lot of accuracy 
is lost at small y with this procedure reasonable accuracy can be achieved at larger values of the 
frequency parameter y. However, high accuracy requires a precise and smooth input Ii[f]{a,b,y') 
together with a fine mesh of y'-values. Considering how fast and easy the Filon-Simpson weights can 
be generated this procedure does not offer real advantages and is not recommended. 
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Figure 6: Absolute value of the relative deviation between numerical and exact result if the relation 
(4.9) is used to obtain the integral hifo] from /i[/o]. The latter was calculated by the Filon-Simpson 
routine with 288 points + asymptotic contribution as function of y' in steps of Ay' = 0.5 and then 
numerically integrated over y' using either the trapezoidal rule or Simpson's rule. For comparison the 
result from the direct evaluation of /2[/o] by the Filon-Simpson rule with 288 integration points + 
asymptotic contribution (as in Fig. 4) is also shown. 



5 Summary 

Relatively simple and straightforward quadrature rules of Filon-Simpson form have been presented 
which are applicable for numerical integration of oscillatory integrals of the type (1.2, 1.3). They 
employ equidistant integration points (including the endpoints) and weights which have to be calcu- 
lated anew for each value of the frequency parameter y. The choice of equidistant points allows easy 
construction of extended Filon-Simpson quadrature rules so that the accuracy of the result can be 
simply assessed by increasing the number of subdivisions. Inevitably the Filon-Simpson weights are 
more involved than the ones of standard quadrature rules as they are given in terms of sine and cosine 
integrals and elementary functions. However, the price for an accurate evaluation of these weights 
is modest and worthwhile as the Filon-Simpson quadrature rules not only reduce to the ordinary 
Simpson rule for y = but also give the leading and subleading terms for ?/ — >^ oo provided smooth 
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functions are integrated and the spacing of integration points is fine enough. Although a rigorous 
error estimate has yet to be given, numerical tests have shown that in this regime they do far better 
than standard quadrature rules. This makes the Filon-Simpson method attractive for all applications 
where these types of rapidly-oscillating integrals have to be evaluated. 



Acknowledgement: I would like to thank Prof. A. Iserles for a very kind correspondence. 



Appendix 



Al Exact asymptotic behaviour of the integrals 



Here we derive the asymptotic expansion of the oscillatory integrals beyond the leading order which 
was only considered in the main text. We assume that the function /(x) is analytic at x = and that 
the upper limit b is finite. 

We start with the integral Ii whose asymptotic expansion is easy to obtain by a subtraction 
followed by an integration by parts so that an additional inverse power of y is generated: 



,[/](a = 0,6,2/) = /' dx /(O) + - f'dx ^^""^ ^^"^ sin xy = ^F^'' (by) 

Jo xy y Jo ™ " 



y 



+- 



y 



fix) - /(O) cosxy 
X y 



y Jo V X ) 



cosxy 



• (A.l) 



Repeating the process in the remaining integral it is seen that it is of higher order. By using the 
exphcit expression for FQ^-*(6y) given in eq. (2.11) together with the asymptotic expansion (2.18) we 
thus obtain 



/i[/](a = 0,6,y) 



y^oo 



2y 



y 



cos yb 



+ 



sin by 



(A.2) 



Similarly we obtain for the integral I2 [/] (a = 0, 6, y) after two subtractions 



/2[/](a = 0,6,y) = MFP(by) + L^Fi^\by) + \ f dx g{x) {1 - cos xy) (A.3) 

y y y Jo 



y 



where g{x) = (/(x) — /(O) — x/'(0))/x^ is regular at x = 0. Therefore we may apply the procedure of 
repeated integration by parts to obtain for the last term 



rb 2 2 

/ dx g'(x) (1 — cos xy) = dx gix) T5(x)sinxy 

^0 y Jo 



+ ... . 



(A.4) 



13 



Finally by employing the explicit expressions for F^'^\hy) given in eq. (2.12) together with the 
asymptotic expansions (2.18) we obtain 

ima = 0,M) + ^ [;'(0) h,, + - ?i<^ + . . . . (A.5) 

Formally the last term in eq. (A.5) is of next-to-next-to-leading order but it is needed to verify 
the relation (1.4) to next-to-leading order. Note also the appearance of logarithmic terms in the 
asymptotic expansion of I2 ; the constant C{h) in eq. (A.5) is given by 

CW = (. + ln6)/'(0)-M+/^,/(£WW^^. (A.6) 

J Q X 

Three integration by parts in the last integral bring it into the form 

C{b) = 7/'(0)-4?^ + m+[/'(6)-6r(6)] In 6+ j'dx x\nx f'^x) (A.7) 
Jo 

which will be needed for comparison with the extended Filon-Simpson rule. 

It is clear that nothing prevents the procedure to be extended to arbitrary order but for our 
purposes this is not needed. One may check the asymptotic expansions (A. 2, A.5) by applying them 
to the functions fi{x) used for the numerical tests: in the limit 6 ^ 00 one obtains full agreement 
with the first terms of the asymptotic expansion of the exact results (4.2, 4.3). 

A2 Asymptotic behaviour of the Filon-Simpson rules 

How do the Filon-Simpson quadrature rules behave in the asymptotic limit ? To answer this question 
we just have to plug the asymptotic expansions of the functions F-''\by) into the expressions for the 
weights. After some algebraic work we obtain 



l^wl ' fi = --/o + — 
y2 y2 



-f f2 cos by 
ojo r 



+ o { ^^') , a = 0,b = 2h (A.8) 



where fi = f{ih) and 



Sfo = ^[ -3/(0) + 4f{h) - f{2h) ] ^-^ /'(O) - -h'f"{0) + .... (A.9) 

Thus one obtains the correct subleading term of the asymptotic expansion except that the derivative 
of the function at a; = is replaced by its (forward) finite-difference approximation Sfo- 

Similarly, one finds that for the integral I2 [/] the Filon-Simpson quadrature rule has the asymptotic 
expansion 

2 

4''^ fi = -tt/o + 4 [ 7/o • Iny + Cpsib)] + O f^^\ a = 0, b = 2h (A.IO) 

i=o y y \ y / 

where the constant is given by 

CFs{b) = (7 + ln6) 7fo-j + ls^fi. (A.ll) 
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Here 

S'h = ^ [ /(O) - 2/(^) + fi^h) ] r (0) + hf"'{0) + ... = f"{h) + {h^) (A.12) 

is a finite-difference approximation to the second derivative of the function f{x) at x = (or better 
sX X = h). Comparing with the exact next-to- leading term in the asymptotic expansion of l2[f] we 
see that again finite differences are subsituted for derivatives and that the last integral in eq. (A. 6) is 
replaced by 



JO 



dx 



f{x)-f{0)-xf'{0) 



Jo 



dx 



p"io) + ... 



(0) + . . . 



(A.13) 



which is valid for regular functions and increments b = 2h which are small enough. 



A3 Asymptotic behaviour of the extended Filon-Simpson rules 

Let us now investigate how the extended Filon-Simpson rules behave in the limit y ^ oo. To do that 
we need the asymptotic behaviour of the simple rules in intervals with non-zero lower limit. Using 
the quadrature rules in the form 



i=0 



fo-aSfo + ^S^fi 



4^' + 



we obtain for the Filon-Simpson quadrature of Ii [/] 

2 



(1). y- 



1=0 



1 

y2 



fo /2 , 

— COS ay — — cos by 
a b 



a^0,b = a + 2h 



(A.14) 



(A.15) 



and therefore 

TV 



i=0 



(1) f y-^oo tt/o J_ 
2y y^ 



cos(2/iy) 
2h 



+- 



1 

y2 
1 

y 



5fo-f2 

f2 

I 

p^cos((2IV-2)M-|^c<»(2«) 



, cos(2/i'u) cos Ahy 

2h ^ ^' 4:h " 



+ ... 



(A.16) 



Here the first line gives the contribution from the first interval [0, 2h] (see eq. (A.8)), the second line 
the one from the second interval [2h, 4h] and so on. It is seen that in the 1/y^ -terms the contributions 
cancel pairwise and only the ones from the first and the last interval survive. Therefore 



N 

E^^'V, 

1=0 



(1). y- 



2y y2 



~Sfo-{^cos2Nh 
Nh 



+ ... 



(A.17) 



which is again the correct asymptotic result (A. 2) except that the derivative of the function at a; = 
is replaced by the finite difference Sfo- 
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The subleading asymptotic terms for the Filon-Simpson quadrature of hif] are more involved 
because of the logarithmic terms. Prom eq. (A. 14) and the asymptotic behaviour of the integrals jj^^ 
one gets after some algebra 



(2) „ y^oo 



1=0 



- - ^ + 2M2/i + (^f2 - bS^h) In 6 
a b ^ ' 

- (5/0 - Ina 



, a /O, 6 = + 2/1. (A.18) 



Here 



5/2 = ^1/0-4/1+3/2] /^-y/f + 



(A.19) 



is the backward finite-difference approximation for the derivative of the function /(x) at the point 
X2 = 6 = a + 2/t. It is obtained from the (forward) form (A. 9) by the exchange a b. 
Together with the result (A. 10) for the first interval we therefore have 



/2i-l 



i=0 ^ ^ K i=2 

+ (Sfo- ~5f2 + 2MV3) ln(2/i) + (I/at - NhS^fN-i) H^h) 

ln(2i/i) I . (A.20) 



N/2-1 
+2h ^ 

i=2 



+ 2in — 



2h 



2h 



At first sight this looks rather complicated unless one recognizes the sums as (extended) trapezoidal 
rules with stepsize 2h for the corresponding integrals (sec, e.g. cq. (25.4.2) in ref. [4]) 



2h 



-J +92 + 94 + ■■■+ 9N-2 + 



j-Nh 

J dx g{x) + O (h^j . (A.21) 



Furthermore 



/2 — fo h-*o f, 
2h ' 

6f2i— 6f2i 

2h 



fo , 'Sfo- 1/2 + 2MV3 = O [h^) , 



<^^/2i+l - i'^hi-l h^O 

2h ' 



fill 

J2i ■ 



(A.22) 



Therefore in the limit h ^ eq. (A.20) becomes 



N 
i=0 



(2) f y-*oo tt/o _2 

y 



fib) 



/'(O) (7 + Iny) + /'(O) -i^+J^dx f"{x) 



+ [ f'{b) - bf"{b) ] Inb + dxxlnx f"{x) 



(A.23) 



which agrees exactly with the subleading term of eq. (A.5) and the form (A.7) of the constant C(6). 
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